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The crystal size distribution in polynuclear growth is numerically studied using a coupled map 
lattice model. The width of the size distribution depends on c/D, where c is the growth rate at 
interface sites and D is the diffusion constant. When c/D is sufficiently small, the width W increases 
linearly with c/D and saturates at large c/D. Monodisperse square and cubic crystals are obtained 
respectively on square and cubic lattices when c/D is sufficiently small for a small kinetic parameter 
b. The linear dependence of W on c/D in a parameter range of small c/D is explained by the 
eigenfunction for the first eigenvalue in a two-dimensional model and a mean-field model. For the 
mean-field model, the slope of the linear dependence is evaluated theoretically. 


I. INTRODUCTION 

Crystal growth has been extensively studied in various research fields such as applied physics, metallurgy, chemical 
engineering, and mineralogy. It is also studied in statistical physics as a typical nonequilibrium phenomenon. Nucle- 
ation, surface kinetics, and diffusion are important factors in crystal growth [lj. Crystal growth is one of the Stefan 
problems for moving interfaces, which is difficult to analyze mathematically. Various simulation methods have been 
proposed to study crystal growth far from equilibrium. The phase-field model is a partial-differential equation in 
which the solid and melt phases are expressed with a continuous variable corresponding to the order parameter 2, 3], 
We proposed a coupled map lattice model for crystal growth M- A similar type of order parameter is defined on a 
lattice, and coupled maps are used for the time evolution. The numerical simulation of the coupled map lattice is much 
faster than that using the phase-field model. One reason is that the width of the interface between the solid and melt 
phases is only one lattice constant in the coupled map lattice model. On the other hand, the order parameter changes 
smoothly at the interface in the phase-field model owing to the diffusion term in the equation for the order parameter. 
We performed numerical simulation of crystal growth of dendrites and diffusion-limited aggregations (DLA’s) using 
coupled map lattices. These complicated fractal patterns can be easily reproduced in the coupled map lattice model, 
because the surface tension effect can be set to zero in the coupled map lattice model. In previous papers, we studied 
crystal growth starting from a single crystal seed. In this paper, we study polynuclear growth in which there are 
many crystal seeds initially, and investigate the crystal size distribution at the final stationary state. 

Crystal size distribution (CSD) is one of important topics in crystallography. Various theories of CSD were proposed. 
Becker and Doring studied the time evolution of the size distribution in the nucleation process [7:]. Randolph and 
Larson discussed stationary CSDs in the research field of chemical engineering [1]. Marsh applied the CSD theory 
in geology and found an exponential size distribution of mineral crystals deposited from magma (t|. In their model, 
the exponential distribution is derived by the balance of the outflow, the inflow, and the growth of crystals. In the 
synthesis of micro- and nanocrystals, the final size distribution after the crystal growth is important. Various methods 
of characterization and control of the crystal size distribution have been studied. In particular, the control of the 
crystal size distribution have been investigated in colloidal sciences pd], El- The size distributions similar to the 
gamma distribution or the log normal distribution are often observed in experiments. The size distribution P(S) 
is not symmetric around the average value but has a longer tail for large S. It is practically important to produce 
almost monodisperse crystals, because the monodisperse crystals are useful for various applications such as catalysts, 
magnetic fluids, and drug delivery systems. Colloidal crystals can be produced from monodisperse microcolloidal 
particles, which can be applied to various systems such as photonic crystals fl2| . 

The crystal size distribution depends on many factors such as nucleation processes, surface conditions of crystals, 
and various experimental conditions. The CSD problem is still an open problem. For the synthesis of monodisperse 
crystals, it is important that the time scales of nucleation and growth processes are separated, and the size distribution 
of the crystal seeds is almost mono-disperse at the initial stage of the growth process. In this paper, we consider a 
problem that the size distribution of the crystal seeds is almost monodisperse at the initial stage of the crystal growth, 
but the initial positions of the crystal seeds are randomly distributed. The final size distribution changes owing to the 
competition of the diffusion and surface kinetics. Our model might be applied to a two-dimensional crystal growth 
on a substrate in that the positions of crystals are fixed in time. However, the purpose of this study is to understand 
qualitatively a condition wherein the final crystal size distribution becomes almost monodisperse. 
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FIG. 1: Crystal sites at (a) D = 0.01 and c = 1, and (b) D = 0.2 and c = 0.01 in the case that the initial sizes of crystals are 
all 1. Crystal sites at (c) D = 0.01 and c = 1, and (d) D = 0.2 and c = 0.01 in the case that the initial sizes of crystals are 1 
for i < 70 and 9 for i > 70. 


II. TWO-DIMENSIONAL COUPLED MAP LATTICE MODELS FOR POLYNUCLEAR GROWTH 

We first study a two-dimensional coupled map lattice model on a square lattice. There are two variables in our 
coupled map lattice for the solution growth, namely, the order parameter x n (i,j) and the concentration u n (i,j), where 
(i,j) and n denote respectively a lattice site and a discrete time. The order parameter x n (i,j) is 1 at crystal sites, 
0 at solution sites, and between 0 and 1 at interface sites. The time evolution has two stages, namely, diffusion and 
surface kinetics. The diffusion is expressed as 

) = u n (i,j) + D{u n (i + 1 ,j) + u n (i - 1 ,j) + u n (i,j + 1) + u n (ij - 1) - 4 u n (i,j)}, (1) 

where D is a diffusion constant. The surface kinetics is expressed as 

Xn+i(i,j) = x n (i,j)+cbu' n (i,j), 

Mi (i,j) = u' n (i,j) - cbu' n (i,j), (2) 

where c is a parameter for the growth rate. The time evolution Eq. (2) is applied only at interface sites satisfying 
0 < x n (i,j) < 1. If x n (i,j ) goes over a critical value 1, interface sites change into crystal sites and the neighboring 
solute sites are assigned to interface sites. By repeating Eqs. (1) and (2), polynuclear crystal growth occurs. The 
parameter b is 1 at nonflat sites and takes a small value bo at flat sites. The flatness of an interface site is determined 
by counting the number of crystal sites in the four nearest-neighbor sites of the interface site. If the number is 1, & 
is set to be &o, and if the number is larger than 2, b is set to be 1. For sufficiently small bo, square crystals appear, 
because kink sites and concave sites are quickly occupied, but it takes a rather long time for flat sites to grow. The 
parameter b represents an effect of the surface tension on the surface kinetics. In the case of &o = 1, the surface 
tension effect disappears. In Eq. (2), the increase in the order parameter x is equal to the absorption of u at interface 
sites, that is, the conservation law of x + it is satisfied. In our model of Eqs. (1) and (2), only the deposition process 
from solution to crystal is assumed, that is, the inverse process of dissolution is not taken into consideration. 

We have performed numerical simulation of the two-dimensional coupled map lattices. Initially, 10 crystal seeds of 
size 1 are randomly set in a square of L x L = 100 x 100. The initial value of u(i,j) is uniform and set to be uq = 0.1, 
and the parameter b is assumed to be 0.005. Periodic boundary conditions are imposed. After a long-time iteration 
of Eqs. (1) and (2), u becomes zero owing to the diffusion and the absorption of u at interface sites, and the crystal 
growth process is completed. Figure 1(a) shows crystal sites at D — 0.01 and c = 1 at the final time. The crystal 
size is widely distributed. The crystal sizes are small in a region where initial seeds are densely distributed. In our 
model, the randomness originates only from the initial configuration of crystal seeds. Figure 1(b) shows crystal sites 
at D = 0.2 and c = 0.01. Crystals take a clear square shape, and the sizes are almost uniform. As D is large and c is 
small, the concentration u tends to be uniform and the crystal sizes become homogeneous even for the random initial 
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FIG. 2: Size distributions at (a) D = 0.01 and c = 1, and (b) D = 0.2 and c = 0.01. (c) Relationship between W and c for 
c/D = 5 (rhombi) and c/D = 2.5 (pluses), (d) Relationship between W and c/D. 



FIG. 3: Size distributions at (a) D = 0.025 and c = 0.64, and (b) D = 0.2 and c = 0.00125 for Eqs. (1),(4), and (5). (c) 
Relationship between W and c/D. 


configuration of seeds. To show an effect of the initial seed size distribution, the sizes of four crystal seeds at i > 70 
are set to be 3 x 3 = 9, and the sizes of the other seeds are set to be 1. Figures 1(c) and 1(d) show crystal sites at (c) 
D = 0.01 and c = 1 and (d) D = 0.2 and c = 0.01 at a final time. At D = 0.2 and c = 0.01, the crystal sizes for the 
six seeds of size 1 are the same as those in Fig. 1(b), and the crystal sizes for the four larger seeds of size 9 are larger 
by ll 2 — 9 2 = 40 than those in Fig. 1(b). At D = 0.01 and c = 1, the crystal sizes change for all crystals from the 
values in Fig. 1(a), but the change in size from Fig. 1(a) depends on the position. The final crystal size depends on 
the initial seeds, but the dependence of the final size distribution on the initial size distribution is not so simple. We 
investigate only the case that the initial sizes are all 1 hereafter. 

We have calculated the CSD using a larger square lattice of 400 x 400. The size S is defined as the sum of x n {i,j) 
in crystal sites and interface sites around each crystal seed. At crystal sites, £„(*,/) takes a value of 1. Then, S 
represents the area of a two-dimensional crystal. Figures 2(a) and 2(b) show the size distribution at (a) D = 0.01 
and c = 1, and (b) D = 0.2 and c = 0.01. (There are a few crystals that collide with neighboring crystals, but the 
sizes of those crystals are not considered in the size distribution. The collision probability decreases if the system size 
increases for a fixed number of crystal seeds. In the production of colloidal particles, the electric double layer around 
colloidal particles prevents the collision.) The size distribution is rather wide at D = 0.01 and c = 1. The distribution 
is asymmetric and has a longer tail for large S. The size distribution is very narrow at D = 0.2 and c = 0.01, which 
implies that crystals are almost monodisperse. 

The average of crystal size is estimated as uqL 2 /N , because the total sum of S is equal to L 2 u o owing to the 
conservation law of u + x. The inhomogeneity of the crystal size is characterized by the width W or the standard 
deviation of the size distribution. Figure 2(c) shows the widths W for four different c’s at the same c/D = 2.5 and 
5. The width W depends strongly on c/D but hardly depends on c under the condition of c/D =const. If c and D 
are small and the one-step increments of x and u in Eqs. (1) and (2) are sufficiently small, Eqs. (1) and (2) can be 
approximated by differential equations. It can be shown that the final stationary states depend only on the ratio c/D 
owing to a scale transformation of time if the differential approximation is good. Figure 2(d) shows a relationship 
between W and the ratio c/D in a logarithmic plot. The diffusion constant D is changed as 0.01,0.025,0.05,0.1 
and 0.2 for c = 1, and then c is changed as 1/2" (n = 1,2, ••• ,6) for D = 0.2. The width increases linearly as 
W ~ 10.6(c/H) in a range of small c/D and tends to saturate at large c/D. That is, monodisperse crystals are 
obtained in the limit of c/D = 0, and the dispersion is proportional to c/D when c/D is sufficiently small. 

In the coupled map lattice model of Eqs. (1) and (2), there is a possibility that neighboring crystals become closer 
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to each other and collide with each other as they grow. To neglect the effect of the collision, we can propose a simpler 
coupled map lattice model. Equation (2) is rewritten as 

x n +i (i,j) = x n (i,j)+du' n (i,j)/r n {i,j)x2irr n (i,j), 

u n +i(i,j) = u n+1 / 2 (i,j) ~ c'u' n (i,j)/r n (i,j) x 2nr n (i,j), (3) 

where c' = cb/(2 tt), and x n {i,j) and r n (i,j) = ■s/x n {i,j)/'x are interpreted respectively as the area and radius of 
circular crystals. Here, the growth velocity of crystals of radius r n (i, j) is assumed to be inversely proportional to the 
radius r n (i,j) as v r = d/r. This type of relation is satisfied under the condition of diffusion-limited growth. Equation 
(3) is expressed as 


x n +i(i,j) = x n (i,j)+cu' n (i,j), (4) 

u n +i(i,j) = u' n (i,j)-cu' n (i,j), (5) 

where 27 tc' is rewritten as c. The form of Eqs. (4) and (5) is almost equal to that of Eq. (2); however, the time 
evolution of Eqs. (4) and (5) is applied only at the positions of crystal seeds in this model. That is, in contrast 
to Eq. (2), the variable x increases indefinitely even if x becomes larger than 1. Growth patterns such as square 
patterns do not appear, because x grows only at the seed points. The magnitude of x is interpreted as crystal size. 
In other words, we consider a situation that the radii of crystals are sufficiently small in contrast to the distances 
among different crystals, and the growth patterns are invisible. In the following numerical simulation, the number N 
of initial seeds is set to be 180, and the system size is Lx L = 400 x 400. The crystal growth stops for sufficiently large 

n, when u n (i,j) decreases to 0. Figure 3(a) shows the crystal size distribution at D = 0.025 and c = 0.64. Figure 

3(b) shows the size distribution at D = 0.2 and c = 0.00125. When c is small and D is large, the size distribution 
becomes narrow. The rhombi in Fig. 3(c) show the relationship between the width W of the size distribution and 
the ratio c/D, where c is changed as 5.12/2" (n = 1, 2, • • • , 12) for D = 0.2. The dashed line denotes 74 (c/D). The 
width W increases from 0 linearly in a parameter region of small c/D and tends to saturate at large c/D also in this 
simplified model. 

Because Eqs. (1) and (5) are linear equations, an initial value problem for u can be solved using the eigenvalues and 
eigenfunctions as u n (i,j) = Y/k a,kek(i,j)( A*,)", where Xk and ek(i,j ) are eigenvalues and eigenfunctions, and a*, is the 
expansion coefficient. The first eigenvalue Ai and the corresponding eigenfunction Ci(i, j) can be numerically evaluated 
by a long-time iteration of the same equation Eqs. (1) and (5) and the normalization u n +i(i, j) —»• u„+i(i, j)/X\. When 
c/D is sufficiently small, the first eigenvalue is close to 1 and the first component of a\e\(i, j)(Xi) n decays slowly or 
the other components a kCk(i, j)(^k) n decay quickly. Under this condition, the growth rate of each crystal is 

expected to be proportional to the first eigenfunction e\ (i,j) at the crystal-seed point, because the concentration u 
decays as e±(i, j)A" and x(i,j ) increases to compensate for the decrease in u. The crystal-size distribution can be 
evaluated approximately using numerically estimated e\(i,j) by assuming that the crystal size is proportional to the 
eigenfunction e\ (i,j) and using the fact that the average crystal size is equal to L 2 uq/N. The pluses in Fig. 3(c) show 
the width of the size distribution evaluated from ei(i, j). The width is close to the results denoted by rhombi obtained 
by direct numerical simulation of Eqs. (1), (4), and (5) for small values of c/D. At large c/D, the contribution from 
the other eigenfunctions cannot be neglected, and some deviation is observed. 


III. VORONOI TESSELLATION AND A MEAN-FIELD MODEL OF POLYNUCLEAR GROWTH 


The Voronoi tesselation is a method of space partitioning. In the construction of the Voronoi tessellation, N central 
points are randomly distributed, and the whole plane is partitioned into territories of each point. The boundary line 
between two territories is determined by the perpendicular bisector of the two central points. Each Voronoi cell is a 
set of points surrounded by perpendicular bisectors. The Voronoi tesselation is applied to various areas such as grains 
of crystals and territories of animals. The size distribution of Voronoi cells is approximately given by the gamma 
distribution: dmi 


P(S V ) oc S/e 


—aS v 


( 6 ) 


where S v is the area of each Voronoi cell, a is a parameter, and 7 ~ 2.5. 

In our model, the crystal grows fast by absorbing the surrounding solute when c/D is sufficiently large. The 
absorption range or the territory of each crystal seed might be approximated by the Voronoi cell. We interpret that 
the central point in each Voronoi cell corresponds to a seed point of a crystal. First, we consider a case that all 
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FIG. 4: (a) Size distribution determined by the Voronoi tessellation where a is given by Eq. (9). (b) Size distribution determined 
using Eqs. (1), (10), and (11) at D = 0.01 and c = 1. (c) Size distribution determined using the mean-field model at D = 1 
and c = 0.02. (d) Relationship between the width W and c/D. 


the solutes in each Voronoi cell are absorbed into the central seed of a crystal. In that case, the crystal size S is 
proportional to the area S v of the Voronoi cell. In that case, the crystal-size distribution is expressed as 

P(S) oc S^e- aS . (7) 


The parameter a is determined from the average value of S. In our model, the average size of crystals is given by 
uqL 2 /N , because of the conservation of the total mass of solute. At 7 = 2.5, the average value can be calculated as 


= f 0 °° S^e-“ s dS = i r( 7 +l) 
' ; f™S'y e -« s dS a r( 7 ) 

from the definition of the gamma function. The parameter a is given by 


7 N 
2 u 0 L 2 


0.0393 
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2a’ 


( 8 ) 

( 9 ) 


for N = 180, Mo = 0.1, and L = 400. Figure 4(a) shows the distribution of Eq. (7) for N = 180, Mo = 0.1, and L = 400. 
The size distributions in Figs. 2(a) and 3(a) are qualitatively close to the size distribution shown in Fig. 4(a), in that 
there is a longer tail in a region of S > (S ). The average value in the size distribution in Fig. 3(a) is close to that 
in Fig. 4(a) because the number of initial crystal seeds is the same; however, the width of the size distribution in 
Fig. 3(a) is smaller than that of the size distribution of the Voronoi tessellation. 

If c/D is not so large, the effect of diffusion needs to be taken into consideration. To understand qualitatively 
the dependence of the width of the size distribution on c/D using an even simpler model, we propose a mean-field 
type model for crystal growth, keeping the Voronoi tessellation in mind. The model is coupled differential equations 
expressed as 


dxi 

Ht = CUi ’ 

^ = -£-{-cUi +D(u-Ui)}, 

where Xi denotes the crystal size in the zth Voronoi cell, iq is the concentration of the solute in the zth Voronoi cell, 
u = (1 /N)J2m, and S v * denotes the size of the zth Voronoi cell. The crystal is assumed to grow only at the seed 
point located at the center of the zth Voronoi cell with the growth rate CUi , which is similar to the previous model 
using Eqs. (1), (4), and (5). The concentration inside each Voronoi cell is assumed to be uniform and denoted as 
Ui. The concentrations iq in different Voronoi cells are different but becomes uniform owing to the mean-field-type 
diffusion term D{u — u/) in Eq. (11). The total mass Y^ii x i + u i^vi) is conserved in the time evolution of Eqs. (10) 
and (11). In our numerical simulation, S vl is determined by the Voronoi tessellation of 180 seeds in a 400 x 400 square 
lattice. The crystal size Si is defined as Xi at the final stationary state of Eqs. (10) and (11). Figure 4(b) is the size 
distribution of Si at the final time at c = 1 and D = 0.01. A rather wide size distribution close to the distribution 
shown in Fig. 4(a) is obtained at this parameter set satisfying c/D = 100. Figure 4(c) is the size distribution of Si 
at c = 0.02 and D = 1. The size distribution is very narrow, when c/D is sufficiently small. Figure 4(d) shows the 
relationship between the width W of the size distribution and c/D. Similarly to previous models, we find the linear 
dependence of W on c/D in the small range of c/D and the saturation at large c/D. 


( 10 ) 

( 11 ) 
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The concentration tt,; can be expanded using eigenvalues A*, (k = 1,2, • • • , N) and the corresponding eigenfunctions 
ek,i of the linear equation (11) as «,(t) = Ek=i a fc e M e_Afe< where ctfc’s are the expansion parameters. Substitution of 
this expression into Eq. (11) yields 


^ k.i — 




N(c - A k S vl + D )' 

By the summation of Eq. (12) from z = 1 to IV, the eigenvalue A satisfies 

N 

N 


-T 

N ^ 


D 


i=i 


c + D — A S v 


= 1 . 


( 12 ) 


(13) 


This is the IVth-order algebraic equation and there are N solutions corresponding to the N eigenvalues. If c = 0, 
A = 0 is a solution to Eq. (13). For sufficiently small c and small A, Eq. (13) is rewritten as 



A_L 

D N 



1=1 


using the Taylor expansion. The first eigenvalue A is therefore expressed as 



(14) 


where S v = (1 /N) EjLi Svj- If c is sufficiently small, the first eigenvalue Ai is small, and tq can be approximated at 
Ui = aqei ,ie~ Xlt because the other terms decay quickly. Owing to Eq. (12), eq* is proportional to l/(c — AiS,,* + D), 
and Ui(t) is expressed as 


Ui(t) 


_ U *Q , A;< 

C “f D eSvi j Sy 


(15) 


Since u(0) = uo, Uio = uo(c + D — cS V i/S v ) ~ uqD for a sufficiently small c. The final stationary value of Xi is given 
by the time integration of cui(t) as 


Si = 


CUiQ 


—Ay t 


S v uqD 


C ~ I - D cSy-i j Sy 


C ~ I - D cSy-i j Sy 


SyUQ ( 1 + 


C Syj So 


D 


(16) 


The average value of S', is S v uq and the width W of the size distribution of Si is given by 

W = A ss vUo , (17) 

where SS V is the root mean square of the distribution of S',„;. The width W is evaluated as W = 47.5 (c/D), since 
SS v Uq is evaluated as y^2/7(L 2 /Nq)uq = 47.5 using a formula of the gamma function. The theoretical result is shown 
by a dashed line in Fig. 4(d). Fairly good agreement is seen in a region of sufficiently small c/D. We have shown 
again that the eigenfunction for the largest eigenvalue determines the width of the size distribution for sufficiently 
small c/D. For sufficiently large c/D, the total solute uoS V i in the zth Voronoi cell is absorbed to the itli seed at the 
center. This is the special situation corresponding to Fig. 4(a). The crystal size is evaluated as S V iUo- Then, the 
width of the size distribution is expected to be W = SSuq = 47.5, which is plotted by the dotted line in Fig. 4(d). 

In this simple model, Xi increases in proportion to Ui and the dynamics of u t is determined only by {rtj}. Therefore, 
if the initial sizes Xi(0 ) are not zero but randomly distributed, the final crystal sizes Xif are calculated as Xif 
- + Xi( 0), where Xifo is the final crystal size in the case of ccj(0) = 0. That is, the final size distribution is 
determined by the convolution of the initial size distribution of x t (0) and the distribution of Xi /o. 


IV. THREE-DIMENSIONAL MODEL 

In this section, we study a three-dimensional coupled map lattice, which is more realistic than the two-dimensional 
model. The three-dimensional extension of Eq. (1) on a cubic lattice is expressed as 
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FIG. 5: (a) 3D plot of the crystalized site at D = 0.15 and c = 0.1 for Eqs. (18) and (19). (b) Size distribution at D = 0.01875 
and c = 1. (c) Size distribution at D = 0.15 and c = 0.125. (d) Relationship between the width W and c/D for c = 1. 


u' n (i,j,k) = u n (i,j,k) + D{u n (i + 1 ,j,k) + u n (i - 1 ,j,k) + u n (i,j + 1 ,k) +u n (i,j - 1 ,k) 

+u n (i,j,k + 1) + u n (i,j,k - 1) - 6 u n (i,j,k)}. (18) 

The three-dimensional extension of Eq. (2) is expressed as 

Xn+t{i,j,k) = x n (i,j, k) + cbu' n (i,j, k), 

u n +\{i,j,k) = u' n (i,j, k) — cbu' n (i,j, k), (19) 

where (i,j,k) denotes a lattice point in a cubic lattice, and the parameter b is 1 for nonflat sites and takes a small 
value bo = 0.005 for flat sites. Here, the flat interface site is a site that has only one crystal site and five solution sites 
as its nearest-neighbor sites. The system size is set to be 60 x 60 x 60. Figure 5(a) shows a 3D plot of the crystal site 
obtained in the numerical simulation of the three-dimensional coupled map lattice at D = 0.15 and c = 0.1. Initially, 
11 crystal seed points are randomly distributed. Clear cubic crystals are created, and the sizes of crystals are almost 
uniform. The size distributions are calculated for 8 samples of initial 100 seeds of crystals. Figure 5(b) shows the size 
distribution at D = 0.01875 and c = 1. Figure 5(c) shows the size distribution at D = 0.15 and c = 0.125. Similarly to 
the two-dimensional models, the size distribution becomes wider as c/D increases. Figure 5(d) shows the relationship 
between the width W of the size distribution and c/D. Here, D is changed as D = 0.3/2" (n = 1,2, ■ • • ,5) for c = 1, 
and then c is changed as 1/2" (n = 1, 2, • • • ,5) for D = 0.15. The width increases linearly with c/D, which is similar 
to the results found in the two-dimensional models. 


V. SUMMARY 

We have numerically studied the crystal size distribution in polynuclear growth using coupled map lattice models. 
The size distribution of crystals originates from the initial random configuration of crystal seeds, and the dispersion is 
reduced by the effect of the diffusion. We have found that the width of the size distribution is approximately determined 
by the parameter c/D, where c is the growth rate at the interface sites and D is the diffusion constant. When c/D 
is sufficiently small, the width W increases linearly with c/D and tends to saturate at large c/D. The numerical 
simulations were performed on a square lattice and a cubic lattice. Monodisperse square and cubic crystals are obtained 
when c/D are sufficiently small and the kinetic parameter b is sufficiently small. The research and development of 
the manufacturing process of monodisperse crystals, especially monodisperse nanocrystals, are important for various 
applications. The linear dependence of W on c/D in the parameter range of small c/D is explained by the eigenfunction 
for the first eigenvalue in a simplified two-dimensional model and an even more simplified mean-field-type model. For 
the mean-field model using the Voronoi tessellation, the slope of the linear dependence has been evaluated theoretically. 
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